Locally adaptive tree-based thresholding using the treethresh package in R

Suppose we have, after possible rescaling to obtain unit variance, observed a sequence X = (Xi)i∈I satisfying Xi = μi + i, for i ∈ I, where μ = (μi)i∈I is a possibly sparse signal (i.e. some/most of the μi are believed to be zero), the i are independentN(0, 1) noise, and I is a possibly multidimensional index domain. Being a generalisation of the EbayesThresh method, the TreeThresh method is based on assuming a mixture between a point mass at zero (denoted δ0) and a signal with density γ(·) as prior distribution for the μi: fprior(μi) = (1− wi)δ0 + wiγ(μi) In contrast to the EbayesThresh method the mixing weights wi depends on the index i, i.e. the underlying signal can be heterogeneous (in the sense of not being everywhere equally sparse). We assume there is a partition of the index space I = P1 ∪ . . . ∪ Pp , Pk ∩ Pl = ∅, such that the weights within each region P are (almost) constant. The treethresh software uses a double exponential distribution1 with fixed scale parameter a (set to 0.5 by default) as γ(·), which is also the default setting used in the EbayesThresh package. This yields